## Preliminaries: Load packages and set working directory ----
if(!require(pacman)) {install.packages("pacman");require("pacman")} # Copied from the class script for polynomial regression
p_load("BayesFactor", "sjmisc", "sjPlot", "effects", "berryFunctions", "rstudioapi", "car", "egg", "tidybayes", "devtools", "rstanarm", "bayestestR", "emmeans", "bayesplot", 
       "ppcor", "lubridate", "see", "XNomial", "tidyverse")

current_path <- getActiveDocumentContext()$path    ## gets path for current R script
setwd(dirname(current_path))                       ## sets dir to R script path

## Load cleaned and shaped data ----

cu.2cups.long <- readRDS(file = "two_cups_data.RDS") # 2 slides data
mc.cu.df <- readRDS(file = "three_and_four_cups_data.RDS") # 3 and 4 options data, both Mody and Carey (2016) and the current experiment

cu.2cups.bayesian.mem <- readRDS(file = "bayesian_2cups_model.RDS") # Load saved 2-options model to avoid lengthy fitting process
mc.cu.bayesian.mem <- readRDS(file = "mc.cu.bayesian.mem.RDS") # Load saved 3- and 4-options model to avoid lengthy refitting

g.df <- readRDS("g.df.RDS") # Load cleaned and shaped data from Grigoroglou et al.

## Main analysis of 2-slides data ----

# Descriptive proportions correct by participant
prop.correct.2cups <- cu.2cups.long %>%
  group_by(id, age.group) %>%
  summarize(prop.correct = mean(response, na.rm = TRUE))

# Descriptive proportion correct by age group
(mean.correct.2cups.by.age <- prop.correct.2cups %>%
    group_by(age.group) %>%
    summarize(prop.correct = mean(prop.correct)))

# Extract the posteriors
cu.2cups.bayesian.emm <- emmeans(cu.2cups.bayesian.mem, specs = "age.group", type = "response") # Extract estimates
cu.2cups.post <- insight::get_parameters(cu.2cups.bayesian.emm) # Get posterior for each parameter of interest (10,000 samples)

# Prepare data for plotting
cu.2cups.emm.plot <- as.data.frame(cu.2cups.bayesian.emm)
cu.2cups.emm.plot$age.group <- factor(cu.2cups.emm.plot$age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs"))

cu.2cups.post.long <- cu.2cups.post %>%
  pivot_longer(everything(), names_to = "age.group", values_to = "sample")
# Plot data with 95% credible intervals, along with a small dot representing each participant's proportion correct
ggplot() +
  ggtitle("2-slides task") +
  #geom_violinhalf(data = cu.2cups.post.long, aes(x = age.group, y = sample), fill = "darkgrey", alpha = .5) +
  geom_dotplot(data = prop.correct.2cups, mapping = aes(x = age.group, y = prop.correct),
               binaxis = "y", stackdir = "center", dotsize = .5, 
               binwidth = 1/48, color = "white", fill = "black") +
  geom_errorbar(data = cu.2cups.emm.plot, mapping = aes(x = age.group, ymin = lower.HPD, ymax = upper.HPD), width = .05) +
  geom_point(data = cu.2cups.emm.plot, mapping = aes(x = age.group, y = prob),
             stat = "identity", alpha = .9, size = 3, shape = 23, fill = "white") +
  scale_shape_manual(values = c(1, 2)) +
  ylim(0,1.05) +
  theme_bw() +
  theme(text = element_text(size=24, family = "serif"), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), plot.background = element_blank(), 
        panel.background = element_blank(), legend.title = element_blank(), strip.background = element_blank()) +
  ylab("Median probability choosing target, 95% CIs \n Small dots: proportion choice of target") +
  xlab("") 

# Calculate the proportion of each posterior that is > .50
length(cu.2cups.post[,1][cu.2cups.post[,1] > .50]) / length(cu.2cups.post[,1]) # 99.23% > .50
length(cu.2cups.post[,2][cu.2cups.post[,2] > .50]) / length(cu.2cups.post[,2]) # 100% > .50
length(cu.2cups.post[,3][cu.2cups.post[,3] > .50]) / length(cu.2cups.post[,3]) # 100% > .50

emmeans(cu.2cups.bayesian.mem, specs = "age.group") %>% pairs()

## Main analysis of 3- and 4-options data ----

# Calculate descriptive proportion correct for each participant

prop.correct.df <- mc.cu.df %>%
  group_by(id, age.group, cups.or.channels, trial.type) %>%
  summarize(prop.correct = mean(response, na.rm = TRUE)) %>%
  mutate(trial.type = factor(trial.type, labels = c("3 options", "4 options")),
         cups.or.channels = factor(cups.or.channels, levels = c("cups", "slides")))

prop.correct.df$age.group <- factor(prop.correct.df$age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs")) # For some reason "mutate" would not let me recode these factor labels, don't know why

# Extract the posteriors
mc.cu.bayesian.distributions.by.study <- emmeans(mc.cu.bayesian.mem, specs = "cups.or.channels", by = c("trial.type", "age.group"), type = "response") # Extract estimates
mc.cu.posteriors <- insight::get_parameters(mc.cu.bayesian.distributions.by.study) # Get posterior for each parameter of interest (10,000 samples)

post.probs <- mc.cu.posteriors %>% # Shape posterior samples for easier analysis
  pivot_longer(cols = everything(), names_to = "condition", values_to = "sample") %>%
  mutate(age.group = factor(substr(condition, nchar(condition), nchar(condition)), labels = c("2.5 yrs", "3 yrs", "4 yrs")),
         trial.type = factor(ifelse(grepl("3options", condition), "3 options", "4 options")),
         cups.or.channels = factor(ifelse(grepl("cups", condition), "cups", "slides")))

# Effect of trial type by age group
emmeans(mc.cu.bayesian.mem, specs = "cups.or.channels", by = c("trial.type", "age.group")) %>% pairs()
# Shape data for plotting
mc.cu.bayesian.emm.plot <- as.data.frame(mc.cu.bayesian.distributions.by.study)
mc.cu.bayesian.emm.plot$age.group <- factor(mc.cu.bayesian.emm.plot$age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs"))
mc.cu.bayesian.emm.plot$cups.or.channels <- factor(mc.cu.bayesian.emm.plot$cups.or.channels, levels = c("cups", "slides"))
mc.cu.bayesian.emm.plot$trial.type <- factor(mc.cu.bayesian.emm.plot$trial.type, levels = c("3options", "4options"), labels = c("3 options", "4 options"))

# Plot the data

ggplot() +
  geom_line(data = mc.cu.bayesian.emm.plot, mapping = aes(x = cups.or.channels, y = prob, group = age.group), linetype = "dashed", color = "grey33") +
  geom_point(data = mc.cu.bayesian.emm.plot, mapping = aes(x = cups.or.channels, y = prob, color = cups.or.channels),
             position = position_dodge(width = .7), stat = "identity", size = 4, shape = 18) +
  scale_shape_manual(values = c(1, 2)) +
  geom_errorbar(data = mc.cu.bayesian.emm.plot, aes(x = cups.or.channels, ymin = lower.HPD, ymax = upper.HPD, group = cups.or.channels, width = .5, color = cups.or.channels),
                position = position_dodge(width = .7), width = .05) +
  geom_dotplot(data = prop.correct.df, mapping = aes(x = cups.or.channels, y = prop.correct, color = cups.or.channels, fill = cups.or.channels),
               binaxis = "y", stackdir = "center", dotsize = .4, alpha = .1, binwidth = .08) +
  facet_wrap(~ age.group + trial.type, nrow = 3, labeller = label_wrap_gen(multi_line=FALSE)) +
  ylim(0,1) +
  scale_fill_brewer(palette = "Dark2") +
  scale_color_brewer(palette = "Dark2") + 
  theme_bw() +
  theme(text = element_text(size=24, family = "serif"), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), 
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_blank(), legend.position = "none") +
  labs(title = "3- and 4-options tasks", x = "", y = "Median probability correct, 95% Cred. Intervals") +
  scale_x_discrete(labels = c("cups", "slides")) +
  xlab("") 

# Simplified figure

ggplot() +
  #geom_line(data = mc.cu.bayesian.emm.plot, mapping = aes(x = cups.or.channels, y = prob, group = age.group), linetype = "dashed", color = "grey33") +
  scale_shape_manual(values = c(1, 2)) +
  geom_errorbar(data = mc.cu.bayesian.emm.plot, aes(x = cups.or.channels, ymin = lower.HPD, ymax = upper.HPD, group = cups.or.channels, width = .5),
                position = position_dodge(width = .7), width = .05) +
  geom_dotplot(data = prop.correct.df, mapping = aes(x = cups.or.channels, y = prop.correct),
               binaxis = "y", stackdir = "center", dotsize = .4, color = "white", fill = "black", binwidth = .08) +
  geom_point(data = mc.cu.bayesian.emm.plot, mapping = aes(x = cups.or.channels, y = prob),
             position = position_dodge(width = .7), stat = "identity", size = 4, shape = 23, color = "black", fill = "white") +
  facet_wrap(~ age.group + trial.type, nrow = 3, labeller = label_wrap_gen(multi_line=FALSE)) +
  ylim(0,1) +
  scale_fill_brewer(palette = "Dark2") +
  scale_color_brewer(palette = "Dark2") + 
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), 
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_blank(), legend.position = "none") +
  labs(title = "A. Summary results", x = "Experiment", y = "Median probability choosing target, 95% CI") +
  scale_x_discrete(labels = c("cups", "slides"))

# Calculate posterior medians
post.medians <- post.probs %>%
  group_by(age.group, trial.type, cups.or.channels) %>%
  summarize(post.median = median(sample))

ggplot(post.probs) +
  facet_wrap(~ age.group + trial.type, ncol = 2, labeller = label_wrap_gen(multi_line=FALSE)) +
  geom_density(aes(x = sample, group = cups.or.channels, fill = cups.or.channels), alpha = .5) +
  scale_fill_brewer(palette = "Dark2") +
  scale_color_brewer(palette = "Dark2") + 
  geom_vline(data = post.medians, aes(xintercept = post.median, color = cups.or.channels)) +
  scale_x_continuous(breaks = c(.25, .5, .75)) +
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), #axis.ticks.y=element_blank(), axis.text.y = element_blank(),
        panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), legend.position = c(.87, .935),
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_rect(fill = "white")) +
  labs(title = "Posterior probability target choice", x = "")

slope.posteriors <- pairs(emmeans(mc.cu.bayesian.mem, specs = "cups.or.channels", by = c("trial.type", "age.group"))) %>% # Extract estimates for the slope parameters in logodds
  insight::get_parameters()  # Get all 10000 samples for each parameter
slope.CIs <- as.data.frame(pairs(emmeans(mc.cu.bayesian.mem, specs = "cups.or.channels", by = c("trial.type", "age.group")))) %>%
  mutate(age.group = factor(age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs")),
         trial.type = factor(trial.type, labels = c("3 options", "4 options")))

# Plot
pivot_longer(slope.posteriors, cols = everything(), names_to = "condition", values_to = "sample") %>% # Shape and plot the data
  mutate(age.group = factor(substr(condition, nchar(condition), nchar(condition)), labels = c("2.5 yrs", "3 yrs", "4 yrs")),
         trial.type = factor(ifelse(grepl("3options", condition), "3 options", "4 options"))) %>%
  ggplot() +
  facet_wrap(~ age.group + trial.type, ncol = 2, labeller = label_wrap_gen(multi_line=FALSE)) +
  geom_density(aes(x = sample), fill = "lightgray") +
  xlim(-2.5, 2.5) +
  geom_segment(data = slope.CIs, aes(x = lower.HPD, y = 0, xend = lower.HPD, yend = .35), color = "darkred") +
  geom_segment(data = slope.CIs, aes(x = upper.HPD, y = 0, xend = upper.HPD, yend = .35), color = "darkred") +
  geom_vline(xintercept = -0.1, linetype = "dashed") +
  geom_vline(xintercept = 0.1, linetype = "dashed") +
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), axis.ticks.y=element_blank(), #axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), 
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_blank(), legend.position = "none") +
  labs(title = "B. Difference between studies", x = "ln(odds in slides/odds in cups)", y = "Density")

tidy.slope.posteriors <- mc.cu.bayesian.mem %>%
  emmeans(specs = "cups.or.channels", by = c("trial.type", "age.group")) %>%
  pairs() %>%
  gather_emmeans_draws() %>%
  rename(sample = .value)

tidy.slope.posteriors.within.exp <- mc.cu.bayesian.mem %>%
  emmeans(specs = "trial.type", by = c("cups.or.channels", "age.group")) %>%
  pairs() %>%
  gather_emmeans_draws() %>%
  rename(sample = .value)

tidy.slope.posteriors %>%
  group_by(contrast, trial.type, age.group) %>%
  summarize(median = median(sample),
            ci = ci(sample, method = "HDI"))

tidy.slope.posteriors.within.exp %>%
  group_by(contrast, cups.or.channels, age.group) %>%
  summarize(median = median(sample),
            ci = ci(sample, method = "HDI"))

# For a more detailed quantificational analysis, I want to have a look at how the ROPE of interest stacks up against all possible ROPEs (or a lot of them, anyway)

slope.post.3opts.2 <- tidy.slope.posteriors %>%
  filter(trial.type == "3options", age.group == "2")
slope.post.4opts.2 <- tidy.slope.posteriors %>%
  filter(trial.type == "4options", age.group == "2")

slope.post.3opts.3 <- tidy.slope.posteriors %>%
  filter(trial.type == "3options", age.group == "3")
slope.post.4opts.3 <- tidy.slope.posteriors %>%
  filter(trial.type == "4options", age.group == "3")

slope.post.3opts.4 <- tidy.slope.posteriors %>%
  filter(trial.type == "3options", age.group == "4")
slope.post.4opts.4 <- tidy.slope.posteriors %>%
  filter(trial.type == "4options", age.group == "4")

# This function takes the posterior and calculates the proportion that is between lb and ub (i.e., lower bound and upper bound)
rope.calculator.3opts.2 <- function(lb, ub){
  length(slope.post.3opts.2$sample[slope.post.3opts.2$sample > lb & slope.post.3opts.2$sample < ub]) / nrow(slope.post.3opts.2)
}
rope.calculator.4opts.2 <- function(lb, ub){
  length(slope.post.4opts.2$sample[slope.post.4opts.2$sample > lb & slope.post.4opts.2$sample < ub]) / nrow(slope.post.4opts.2)
}

rope.calculator.3opts.3 <- function(lb, ub){
  length(slope.post.3opts.3$sample[slope.post.3opts.3$sample > lb & slope.post.3opts.3$sample < ub]) / nrow(slope.post.3opts.3)
}
rope.calculator.4opts.3 <- function(lb, ub){
  length(slope.post.4opts.3$sample[slope.post.4opts.3$sample > lb & slope.post.4opts.3$sample < ub]) / nrow(slope.post.4opts.3)
}

rope.calculator.3opts.4 <- function(lb, ub){
  length(slope.post.3opts.4$sample[slope.post.3opts.4$sample > lb & slope.post.3opts.4$sample < ub]) / nrow(slope.post.3opts.4)
}
rope.calculator.4opts.4 <- function(lb, ub){
  length(slope.post.4opts.4$sample[slope.post.4opts.4$sample > lb & slope.post.4opts.4$sample < ub]) / nrow(slope.post.4opts.4)
}

# Set up a grid of hypotheses, so that we can evaluate the probability of each hypothesis
eval.rope.3opts.2 <- tibble(grid = seq(slope.post.3opts.2 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.3opts.2 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                           lower.rope.bound = grid - .1,
                           upper.rope.bound = grid + .1,
                           rope.prob = NA) 
eval.rope.4opts.2 <- tibble(grid = seq(slope.post.4opts.2 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.4opts.2 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                            lower.rope.bound = grid - .1,
                            upper.rope.bound = grid + .1,
                            rope.prob = NA) 

eval.rope.3opts.3 <- tibble(grid = seq(slope.post.3opts.3 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.3opts.3 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                            lower.rope.bound = grid - .1,
                            upper.rope.bound = grid + .1,
                            rope.prob = NA) 
eval.rope.4opts.3 <- tibble(grid = seq(slope.post.4opts.3 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.4opts.3 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                            lower.rope.bound = grid - .1,
                            upper.rope.bound = grid + .1,
                            rope.prob = NA) 

eval.rope.3opts.4 <- tibble(grid = seq(slope.post.3opts.4 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.3opts.4 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                            lower.rope.bound = grid - .1,
                            upper.rope.bound = grid + .1,
                            rope.prob = NA) 
eval.rope.4opts.4 <- tibble(grid = seq(slope.post.4opts.4 %>%
                                         filter(sample < ci(sample)$CI_low) %>%
                                         pull(sample) %>% # This code chunk finds the highest number that is below the 95% CI
                                         max() %>%
                                         round(3),
                                       slope.post.4opts.4 %>%
                                         filter(sample > ci(sample)$CI_high) %>%
                                         pull(sample) %>% # This code chunk finds the lowest number that is above the 95% CI
                                         min() %>%
                                         round(3), 
                                       by = .001), 
                            lower.rope.bound = grid - .1,
                            upper.rope.bound = grid + .1,
                            rope.prob = NA) 

for(i in 1:nrow(eval.rope.3opts.2)){
  eval.rope.3opts.2$rope.prob[i] <- rope.calculator.3opts.2(eval.rope.3opts.2$lower.rope.bound[i], eval.rope.3opts.2$upper.rope.bound[i])
}
for(i in 1:nrow(eval.rope.4opts.2)){
  eval.rope.4opts.2$rope.prob[i] <- rope.calculator.4opts.2(eval.rope.4opts.2$lower.rope.bound[i], eval.rope.4opts.2$upper.rope.bound[i])
}

for(i in 1:nrow(eval.rope.3opts.3)){
  eval.rope.3opts.3$rope.prob[i] <- rope.calculator.3opts.3(eval.rope.3opts.3$lower.rope.bound[i], eval.rope.3opts.3$upper.rope.bound[i])
}
for(i in 1:nrow(eval.rope.4opts.3)){
  eval.rope.4opts.3$rope.prob[i] <- rope.calculator.4opts.3(eval.rope.4opts.3$lower.rope.bound[i], eval.rope.4opts.3$upper.rope.bound[i])
}

for(i in 1:nrow(eval.rope.3opts.4)){
  eval.rope.3opts.4$rope.prob[i] <- rope.calculator.3opts.4(eval.rope.3opts.4$lower.rope.bound[i], eval.rope.3opts.4$upper.rope.bound[i])
}
for(i in 1:nrow(eval.rope.4opts.4)){
  eval.rope.4opts.4$rope.prob[i] <- rope.calculator.4opts.4(eval.rope.4opts.4$lower.rope.bound[i], eval.rope.4opts.4$upper.rope.bound[i])
}

eval.rope.3opts.2[eval.rope.3opts.2$grid == 0.00,] # The probability of the interval around 0 for 2yos, 3 options is .211. 
quantile(eval.rope.3opts.2$rope.prob, probs = seq(.83, .84, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 83% of hypotheses
eval.rope.3opts.2[almost.equal(eval.rope.3opts.2$grid, 0.5),] # The probability of the interval around .5 for 2yos, 3 options is .125.
quantile(eval.rope.3opts.2$rope.prob, probs = seq(.44, .45, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 44% of hypotheses
nrow(slope.post.3opts.2 %>%
       filter(sample > .9, sample < 1.1)) / nrow(slope.post.3opts.2) # Hand calculation for the ROPE around 1: Probabiltiy is .010. Quantile is NA because it's not in the distribution.
nrow(slope.post.3opts.2 %>%
       filter(sample > median(sample) - 0.1, sample < median(sample) + 0.1)) / 10000 # Hand calculation for the median ROPE: Probabiltiy is .2217. 
quantile(eval.rope.3opts.2$rope.prob, probs = seq(.98, .99, length.out = 5)) # Of the hypotheses in the 95% CI, it's the 1.5%
eval.rope.3opts.2 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is 0.121


eval.rope.4opts.2[eval.rope.4opts.2$grid == 0.00,] # The probability of the interval around 0 for 2yos, 4 options is .203. 
quantile(eval.rope.4opts.2$rope.prob, probs = seq(.7, .72, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 70% of hypotheses
eval.rope.4opts.2[almost.equal(eval.rope.4opts.2$grid, 0.5),] # The probability of the interval around .5 for 2yos, 4 options is .163. 
quantile(eval.rope.4opts.2$rope.prob, probs = seq(.52, .55, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 51% of hypotheses
nrow(slope.post.4opts.2 %>%
       filter(sample > .9, sample < 1.1)) / nrow(slope.post.4opts.2) # Hand calculation: Probability is .015. Quantile is NA because it's not in the distribution.
nrow(slope.post.4opts.2 %>%
       filter(sample > median(sample) - 0.1, sample < median(sample) + 0.1)) / 10000 # Hand calculation for the median ROPE: Probabiltiy is .2272. 
quantile(eval.rope.4opts.2$rope.prob, probs = seq(.93, .95, length.out = 5)) # In the 93%ile
eval.rope.4opts.2 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is 0.128 and 0.129


eval.rope.3opts.3[eval.rope.3opts.3$grid == 0.00,] # The probability of the interval around 0 for 3yos, 3 options is .212. 
quantile(eval.rope.3opts.3$rope.prob, probs = seq(.95, .96, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 95% of hypotheses
eval.rope.3opts.3[almost.equal(eval.rope.3opts.3$grid, 0.50),] # The probability of the interval around 0.5 for 3yos, 3 options is .086. 
quantile(eval.rope.3opts.3$rope.prob, probs = seq(.27, .3, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 27% of hypotheses
nrow(slope.post.3opts.3 %>%
       filter(sample > .9, sample < 1.1)) / nrow(slope.post.3opts.3) # Hand calculation: Probability is .005. Quantile is NA because it's not in the distribution.
eval.rope.3opts.3 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is 0.026 and 0.027

eval.rope.4opts.3[eval.rope.4opts.3$grid == 0.00,] # The probability of the interval around 0 for 3yos, 4 options is .234. 
quantile(eval.rope.4opts.3$rope.prob, probs = seq(.975, .98, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 97% of hypotheses
eval.rope.4opts.3[almost.equal(eval.rope.4opts.3$grid, 0.50),] # The probability of the interval around 0.5 for 3yos, 3 options is .093. 
quantile(eval.rope.4opts.3$rope.prob, probs = seq(.295, .3, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 29% of hypotheses
nrow(slope.post.4opts.3 %>%
       filter(sample > .9, sample < 1.1)) / nrow(slope.post.4opts.3) # Hand calculation: Probability is .006. Quantile is NA because it's not in the distribution.
eval.rope.4opts.3 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is -0.011

eval.rope.3opts.4[eval.rope.3opts.4$grid > -0.1,] # The probability of the interval around 0 for 4yos, 3 options isn't even on the map
nrow(slope.post.3opts.4 %>%
       filter(sample > -.1, sample < .1)) / nrow(slope.post.3opts.4) # Hand calculation: Probabiltiy is .007. Quantile is NA because it's not in the distribution.
eval.rope.3opts.4 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is -0.969


eval.rope.4opts.4[eval.rope.4opts.4$grid == 0.00,] # The probability of the interval around 0 for 4yos, 4 options is .136. 
quantile(eval.rope.4opts.4$rope.prob, probs = seq(.48, .49, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 48% of hypotheses
eval.rope.4opts.4[almost.equal(eval.rope.4opts.4$grid, 0.50),] # The probability of the interval around 0.5 for 4yos, 4 options is .207. 
quantile(eval.rope.4opts.4$rope.prob, probs = seq(.75, .79, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 78% of hypotheses
eval.rope.4opts.4[almost.equal(eval.rope.4opts.4$grid, 1),] # The probability of the interval around 1 for 4yos, 4 options is .043. 
quantile(eval.rope.4opts.4$rope.prob, probs = seq(.05, .075, length.out = 5)) # Of the hypotheses in the 95% CI, it's better than 4% of hypotheses
eval.rope.4opts.4 %>% filter(rope.prob == max(rope.prob)) # Most likely hypothesis is 0.368


ggplot(eval.rope.4opts.2, aes(x = grid, y = rope.prob)) +
  geom_bar(stat = "identity", color = "lightgrey", fill = "lightgrey") +
  ylim(0, .24) +
  theme_bw()

# Exploration of weird 4-year-old result
emmeans(mc.cu.bayesian.mem, specs = "age.group", by = c("trial.type", "cups.or.channels")) %>% pairs(reverse = TRUE) # trial.type = 3options, cups.or.channels = slides, age.group4 - age.group3: logOR -0.49, 95% CI [-1.11, 0.13]
emmeans(mc.cu.bayesian.mem, specs = "trial.type", by = c("age.group", "cups.or.channels")) %>% pairs(reverse = TRUE) # age.group - 4, cups.or.channels = slides, 4options - 3options: logOR 0.97, 95% CI [0.44, 1.47]


## Further analysis: Effect of trial number ----
# set.seed(7492)
 # mc.cu.trial.num.mem <- stan_glmer(response ~ cups.or.channels * trial.type * age.group * trial.number + (1|id),
 #                                   family = "binomial", data = mc.cu.df, iter = 3500, warmup = 1000)
 # saveRDS(mc.cu.trial.num.mem, "mc.cu.trial.num.mem.RDS")
mc.cu.trial.num.mem <- readRDS("mc.cu.trial.num.mem.RDS")

mc.cu.trial.num.means <- emtrends(mc.cu.trial.num.mem, specs = ~ cups.or.channels + trial.type + age.group, var = "trial.number", type = "response")
mc.cu.trial.num.post <- gather_emmeans_draws(mc.cu.trial.num.means)

mc.cu.trial.num.medians <- mc.cu.trial.num.post %>%
  group_by(cups.or.channels, age.group, trial.type) %>%
  summarize(median = median(.value))

mc.cu.trial.num.cis <- mc.cu.trial.num.post %>%
  group_by(cups.or.channels, age.group, trial.type) %>%
  summarize(ci = ci(.value, method = "HDI")) %>%
  unnest(ci)

trial.num.plot <- emmip(mc.cu.trial.num.mem, cups.or.channels + trial.type + age.group ~ trial.number, at = list(trial.number = seq(1, 6, by = .5)), CIs = TRUE, plotit= FALSE, type = "response")
all.trial.num.props <- mc.cu.df %>%
  group_by(age.group, trial.type, trial.number) %>%
  summarize(prop.correct = mean(response))

trial.num.props <- mc.cu.df %>%
  group_by(age.group, cups.or.channels, trial.type, trial.number) %>%
  summarize(prop.correct = mean(response, na.rm = TRUE)) %>%
  mutate(age.group = factor(age.group, levels = c("2", "3", "4"), labels = c("2.5 yrs", "3 yrs", "4 yrs")))

ggplot(trial.num.plot %>% filter(cups.or.channels == "slides") %>% 
         mutate(age.group = factor(age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs"))), 
       aes(x = trial.number, y = yvar)) +
  facet_wrap(~age.group + trial.type, nrow = 3) +
  geom_smooth() +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = .25, color = "white") +
  geom_point(data = trial.num.props %>% filter(cups.or.channels == "slides"), aes(y = prop.correct)) +
  ylim(0,1) +
  theme_bw() +
  theme(panel.grid.minor = element_blank()) +
  labs(x = "Trial number", y = "Probability of choosing target")

# Same analysis, but cutting the cups data and collapsing trial type to check learning across all 12 trials

# set.seed(7492)
# cu.trial.num.mem <- stan_glmer(response ~ age.group * trial.number + (1|id), family = "binomial", 
#                                data = mc.cu.df %>%
#                                  filter(cups.or.channels == "slides") %>%
#                                  mutate(trial.number = ifelse(trial.type == "4options", trial.number + 6, trial.number)), 
#                                iter = 3500, warmup = 1000)
# saveRDS(cu.trial.num.mem, "cu.trial.num.mem.RDS")
cu.trial.num.mem <- readRDS("cu.trial.num.mem.RDS")

cu.trial.num.trends <- emtrends(cu.trial.num.mem, specs = "age.group", var = "trial.number", type = "response", CIs = TRUE) 
cu.trial.num.post <- gather_emmeans_draws(cu.trial.num.trends)
cu.trial.num.trends <- cu.trial.num.trends %>% as_tibble() %>%
  mutate(age.group = factor(age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs")))
cu.trial.num.plot <- emmip(cu.trial.num.mem, age.group ~ trial.number, at = list(trial.number = seq(1, 12, by = .5)), CIs = TRUE, plotit= FALSE, type = "response") %>% as_tibble() %>%
  mutate(age.group = factor(age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs")))

slides.12trial.props <- mc.cu.df %>%
  filter(cups.or.channels == "slides") %>%
  mutate(trial.number = ifelse(trial.type == "4options", trial.number + 6, trial.number),
         age.group = factor(age.group, labels = c("2.5 yrs", "3 yrs", "4 yrs"))) %>%
  group_by(age.group, trial.number) %>%
  summarize(prop.correct = mean(response))

ggplot(cu.trial.num.plot, aes(x = trial.number, y = yvar)) +
  facet_wrap(~age.group) +
  geom_smooth(color = "black", size = .5) +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = .25, color = "white") +
  geom_point(data = slides.12trial.props, aes(y = prop.correct)) +
  geom_text(data = cu.trial.num.trends %>% mutate(trial.number.trend = paste("\U03B2 = ", trial.number.trend %>% round(3))),
            aes(x = 6, y = .95, label = trial.number.trend)) +
  geom_text(data = cu.trial.num.trends %>% mutate(ci = paste0("95% CI = [", lower.HPD %>% round(2), " , ", upper.HPD %>% round(2), "]")),
            aes(x = 6, y = .85, label = ci)) +
  scale_x_continuous(breaks = seq(1,12)) +
  ylim(0,1) +
  theme_bw() +
  theme(panel.grid.minor = element_blank()) +
  labs(x = "Trial number", y = "Prob. of choosing target, 95% CI", title = "Change over trials")

## Descriptives of posterior samples ----

apply(slope.posteriors, 2, mean) 

apply(slope.posteriors, 2, median) 

apply(slope.posteriors, 2, map_estimate)

# Calculate the probability, for each posterior, of a small effect as defined by Chen et al. 2010
small.effect.calculator <- function(x){
  length(x[x < -.52 | x > .52]) / length(x)
}

apply(slope.posteriors, 2, small.effect.calculator)

# Calculate the probability, for each posterior, of a small positive effect as defined by Chen et al. 2010
small.positive.effect.calculator <- function(x){
  length(x[x > .52]) / length(x)
}

apply(slope.posteriors, 2, small.positive.effect.calculator)

# Calculate the probability, for each posterior, of a medium effect as defined by Chen et al. 2010
medium.effect.calculator <- function(x){
  length(x[x < -1.24 | x > 1.24]) / length(x)
}

apply(slope.posteriors, 2, medium.effect.calculator)

# Calculate the probability, for each posterior, of a medium positive effect as defined by Chen et al. 2010
medium.positive.effect.calculator <- function(x){
  length(x[x > 1.24]) / length(x)
}

apply(slope.posteriors, 2, medium.positive.effect.calculator)

# Calculate the proportion of each slope posterior that is negative 
log.odds.are.negative <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x < 0]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.negative) # Apply the function over the columns in the dataframe

# Calculate the proportion of each slope posterior that is positive
log.odds.are.positive <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x > 0]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.positive) # Apply the function over the columns in the dataframe

# Calculate the proportion of each slope posterior that is > .25
log.odds.are.more.than.25 <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x > 0.25]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.more.than.25) # Apply the function over the columns in the dataframe

# Calculate the proportion of each slope posterior that is > .5
log.odds.are.more.than.5 <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x > 0.5]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.more.than.5) # Apply the function over the columns in the dataframe

# Calculate the proportion of each slope posterior that is > .75
log.odds.are.more.than.75 <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x > 0.75]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.more.than.75) # Apply the function over the columns in the dataframe

# Calculate the proportion of each slope posterior that is > 1
log.odds.are.more.than.1 <- function(x){ # Define a function that can be applied over each column in the dataframe of posteriors
  length(x[x > 1]) / length(x)
}

apply(slope.posteriors, 2, log.odds.are.more.than.1) # Apply the function over the columns in the dataframe


# Calculate the proportion of each slope posterior that is in the ROPE 
# the rope function from bayestestR does this for us

apply(slope.posteriors, 2, rope) # Apply the function over the columns in the dataframe

# Followup comparisons, since the 4yos came out weird:
# First find the cups-slides contrast for 4yos, 3 slides:
pairs(emmeans(mc.cu.bayesian.mem, specs = "cups.or.channels", by = c("trial.type", "age.group")))

# Next find all the pairwise contrasts of 3 slides to 4 slides
# Note: Signs are backwards here; I want 4opts - 3 opts, but it gives us 3opts - 4opts. Just switch the sign when reporting.
followup.4s.slides <- pairs(emmeans(mc.cu.bayesian.mem, specs = "trial.type", by = c("age.group", "cups.or.channels")))

# Next find the contrasts of 3yos to 4yos
# Note: Signs are backwards here; I want 4opts - 3 opts, but it gives us 3opts - 4opts. Just switch the sign when reporting.
followup.3s.4s <- pairs(emmeans(mc.cu.bayesian.mem, specs = "age.group", by = c("trial.type", "cups.or.channels")))

# Probability of some ROPEs of interest

tidy.main.post <- mc.cu.bayesian.mem %>%
  emmeans(specs = "cups.or.channels", by = c("age.group", "trial.type")) %>%
  pairs() %>%
  gather_emmeans_draws() %>%
  mutate(sample = .value)

#  RsOPE around the median
rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "3options"], range = c(0.01, 0.21))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "3options"], range = c(-0.1, 0.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "3options"], range = c(-1.1, -.9))

rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "4options"], range = c(0.1, 0.3))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "4options"], range = c(-0.07, 0.13))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "4options"], range = c(0.25, 0.45))

# RsOPE around 0
rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "3options"], range = c(-0.1, 0.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "3options"], range = c(-0.1, 0.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "3options"], range = c(-0.1, 0.1))

rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "4options"], range = c(-0.1, 0.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "4options"], range = c(-0.1, 0.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "4options"], range = c(-0.1, 0.1))

# RsOPE around .50 
rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "3options"], range = c(.4, 0.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "3options"], range = c(.4, 0.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "3options"], range = c(.4, 0.6))

rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "4options"], range = c(.4, 0.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "4options"], range = c(.4, 0.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "4options"], range = c(.4, 0.6))

# RsOPE around 1 
rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "3options"], range = c(.9, 1.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "3options"], range = c(.9, 1.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "3options"], range = c(.9, 1.1))

rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "4options"], range = c(.9, 1.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "4options"], range = c(.9, 1.1))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "4options"], range = c(.9, 1.1))

# RsOPE around 1.5
rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "3options"], range = c(1.4, 1.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "3options"], range = c(1.4, 1.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "3options"], range = c(1.4, 1.6))

rope(tidy.main.post$sample[tidy.main.post$age.group == 2 & tidy.main.post$trial.type == "4options"], range = c(1.4, 1.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 3 & tidy.main.post$trial.type == "4options"], range = c(1.4, 1.6))
rope(tidy.main.post$sample[tidy.main.post$age.group == 4 & tidy.main.post$trial.type == "4options"], range = c(1.4, 1.6))


## Extract interaction effects for main analysis ----

mc.cu.bayesian.mem %>%
  emmeans(specs = ~ cups.or.channels * trial.type | age.group) %>%
  contrast(interaction = c("pairwise", "pairwise")) %>%
  insight::get_parameters() %>%
  pivot_longer(everything(), names_to = "contrast", values_to = "sample") %>%
  mutate(age.group = factor(substr(contrast, nchar(contrast), nchar(contrast)), labels = c("2.5 yrs", "3 yrs", "4 yrs"))) ->
interaction.posts.by.age

map_estimate(interaction.posts.by.age$sample[interaction.posts.by.age$age.group == "2.5 yrs"])
map_estimate(interaction.posts.by.age$sample[interaction.posts.by.age$age.group == "3 yrs"])
map_estimate(interaction.posts.by.age$sample[interaction.posts.by.age$age.group == "4 yrs"])

# Find probability of a small negative effect
small.negative.effect.calculator <- function(x) {length(x[x < -.52]) / length(x)}
medium.negative.effect.calculator <- function(x) {length(x[x < -1.24]) / length(x)}



# Some summary stats: Posterior mean and median for each interaction of interest
interaction.posts.by.age %>%
  group_by(age.group) %>%
  summarize(age.mean = mean(sample),
            age.median = median(sample),
            age.ci.low = ci(sample, method = "HDI")[2],
            age.ci.high = ci(sample, method = "HDI")[3],
            prob.small.effect = small.effect.calculator(sample),
            prob.medium.effect = medium.effect.calculator(sample),
            prob.small.negative = small.negative.effect.calculator(sample),
            prob.medium.negative = medium.negative.effect.calculator(sample),
            prob.small.positive = small.positive.effect.calculator(sample),
            prob.medium.positive = medium.positive.effect.calculator(sample))


# Plot these interactions

ggplot(interaction.posts.by.age, aes(x = sample)) +
  facet_wrap(~age.group, nrow = 3) +
  geom_density(fill = "lightgray") +
  geom_vline(aes(xintercept = 0), linetype = "dashed") +
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), axis.ticks.y=element_blank(), axis.text.y = element_blank(),
        panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), legend.position = c(.87, .935),
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_rect(fill = "white")) +
  labs(title = "Experiment:trial type interaction in each age group", x = "")



## Analysis for Study 2 ----


cu5.df <- (read.csv("cu5.csv")) %>% # Load in Eimantas' in-lab 4 year old data
  filter(exclude != "y") %>%
  mutate(age.group = 4) 
cu5.df[cu5.df == "certain "] <- "certain" # Fix a typo

# Participant info

cu5.df %>%
  separate(age, into = c("age.months", "age.days"), sep = ";") %>%
  mutate(approx.age.days = as.numeric(age.months) * 30 + as.numeric(age.days),
         approx.age.years = approx.age.days / 365.25) %>%
  select(approx.age.years) %>%
  summary()

ages <- cu5.df %>%
  separate(age, into = c("age.months", "age.days"), sep = ";") %>%
  mutate(approx.age.days = as.numeric(age.months) * 30 + as.numeric(age.days),
         approx.age.years = approx.age.days / 365.25) %>%
  select(approx.age.years) %>%
  tibble()

sd(ages$approx.age.years)

table(cu5.df$gender)

# Find mean proportion correct

cu5.catching.long <- cu5.df %>%
  select(id, pretest1:pretest6) %>%
  pivot_longer(-id, names_to = "trial", values_to = "choice") %>%
  mutate(response = ifelse(choice == "certain", 1, 0),
         trial.type = "3options",
         experiment = "cu5")

cu5.catching.long %>%
  summarize(prop.correct = mean(response))

# Load Mody & Carey data

mc.data <- read.csv("modycareydata.csv")

mc.4s.choices <- mc.data %>%
  filter(age.group == 4) %>%
  select(id = ID, training1:training3) %>%
  pivot_longer(-id, names_to = "trial", values_to = "response") %>%
  mutate(trial.type = "3options",
         experiment = "mc",
         id = id + 1000)

mc.cu4.cu5.4s.long <- rbind(mc.4s.choices, cu5.catching.long %>% filter(trial.type == "3options") %>% select(-choice)) %>%
  mutate(experiment = factor(experiment, levels = c("cu5", "mc")))

# Model the new catching data
# set.seed(842)
# mc.cu4.cu5.catching.b.mem <- stan_glmer(response ~ experiment + (1|id), family = "binomial", data = mc.cu4.cu5.4s.long, iter = 3500, warmup = 1000)
# saveRDS(mc.cu4.cu5.catching.b.mem, "mc.cu4.cu5.catching.b.mem.RDS")
mc.cu4.cu5.catching.b.mem <- readRDS("mc.cu4.cu5.catching.b.mem.RDS")
mc.cu4.cu5.means <- as.data.frame(emmeans(mc.cu4.cu5.catching.b.mem, specs = ~ experiment, type = "response")) %>%
  mutate(cups.or.channels = "slides",
         trial.type = "3options",
         age.group = 4,
         blocker = 0) %>%
  select(experiment, cups.or.channels, trial.type, age.group, prob, lower.HPD, upper.HPD, blocker)

study2.participant.props <- mc.cu4.cu5.4s.long %>%
  group_by(experiment, id) %>%
  summarize(prop.correct = mean(response)) %>%
  mutate(experiment = factor(experiment, levels = c("mc", "cu5"), labels = c("Mody & Carey", "Study 2")))

ggplot(mc.cu4.cu5.means %>% mutate(experiment = factor(experiment, levels = c("mc", "cu5"), labels = c("Mody & Carey", "Study 2")),
                                   trial.type = "3 options"), 
       aes(x = experiment, y = prob)) +
  geom_dotplot(data = study2.participant.props, aes(x = experiment, y = prop.correct),
               binaxis = "y", stackdir = "center", dotsize = 1,
               binwidth = 1/48, color = "white", fill = "black") +
  geom_errorbar(aes(ymin = lower.HPD, ymax = upper.HPD), width = .05) +
  geom_point(shape = 23, fill = "white") +
  ylim(0,1) +
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), 
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_blank(), legend.position = "none") +  
  labs(title = "A. Summary results", y = "Median probability, 95% CI", x = "Experiment")

cu5.posterior <- insight::get_parameters(emmeans(mc.cu4.cu5.catching.b.mem, specs = ~experiment, type = "response")) %>%
  select(cu5) 
median(cu5.posterior$cu5)
str(ci(cu5.posterior, method = "HDI"))

mc.cu4.cu5.post <- insight::get_parameters(pairs(emmeans(mc.cu4.cu5.catching.b.mem, specs = ~experiment))) %>%
  pivot_longer(everything(), names_to = "contrast", values_to = "sample")

mc.cu4.cu5.post.summary <- mc.cu4.cu5.post %>%
  summarize(median = median(sample),
            lower.ci = ci(sample, method = "HDI")$CI_low,
            upper.ci = ci(sample, method = "HDI")$CI_high)

ggplot(mc.cu4.cu5.post) +
  geom_density(aes(x = sample), fill = "lightgray") +
  xlim(-2.5, 2.5) +
  geom_vline(xintercept = -0.1, linetype = "dashed") +
  geom_vline(xintercept = 0.1, linetype = "dashed") +
  geom_segment(data = mc.cu4.cu5.post.summary, aes(x = lower.ci, y = 0, xend = lower.ci, yend = .35), color = "darkred") +
  geom_segment(data = mc.cu4.cu5.post.summary, aes(x = upper.ci, y = 0, xend = upper.ci, yend = .35), color = "darkred") +
  theme_bw() +
  theme(text = element_text(size=18, family = "serif"), axis.ticks.y=element_blank(), #axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(),
        panel.background = element_blank(), legend.title = element_blank(), 
        strip.background = element_blank(), plot.background = element_blank(), legend.background = element_blank(), legend.position = "none") +
  labs(title = "B. Difference between studies", x = "ln(odds in St. 2/odds in M&C (2016))", y = "Density")

ci(mc.cu4.cu5.post, ci = 0.80, method = "HDI")
str(ci(mc.cu4.cu5.post, ci = 0.71, method = "HDI"))

length(mc.cu4.cu5.post$sample[mc.cu4.cu5.post$sample > 0]) / length(mc.cu4.cu5.post$sample)

rope(mc.cu4.cu5.post, range = c(-0.1, 0.1))
rope(mc.cu4.cu5.post, range = c(-0.57, -0.37))

# For a more detailed quantificational analysis, I want to have a look at how the ROPE of interest stacks up against all possible ROPEs (or a lot of them, anyway)

# This function takes the posterior and calculates the proportion that is between lb and ub (i.e., lower bound and upper bound)
rope.calculator.mc.cu5 <- function(lb, ub){
  length(mc.cu4.cu5.post$sample[mc.cu4.cu5.post$sample > lb & mc.cu4.cu5.post$sample < ub]) / nrow(mc.cu4.cu5.post)
}

# Set up a grid of hypotheses, so that we can evaluate the probability of each hypothesis
eval.rope.mc.cu5 <- tibble(grid = seq(min(mc.cu4.cu5.post$sample) %>% round(2), max(mc.cu4.cu5.post$sample) %>% round(2), by = .001), 
                                lower.rope.bound = grid - .1,
                                upper.rope.bound = grid + .1,
                                rope.prob = NA) 

for(i in 1:nrow(eval.rope.mc.cu5)){
  eval.rope.mc.cu5$rope.prob[i] <- rope.calculator.mc.cu5(eval.rope.mc.cu5$lower.rope.bound[i], eval.rope.mc.cu5$upper.rope.bound[i])
}

eval.rope.mc.cu5[eval.rope.mc.cu5$rope.prob == 0.1183,] # The probability of the interval around 0 is .118. Note that this is not what's reported in the paper, since I used the full posterior here, but it doesn't matter
eval.rope.mc.cu5[eval.rope.mc.cu5$grid == .47,] # Here's the probability of the ROPE around the median
quantile(eval.rope.mc.cu5$rope.prob, probs = seq(.77, .78, length.out = 5)) # Of all hypotheses, it's better than 77% of hypotheses

ggplot(eval.rope.mc.cu5, aes(x = grid, y = rope.prob)) +
  geom_bar(stat = "identity", color = "lightgrey", fill = "lightgrey") +
  geom_vline(xintercept = 0.1, linetype = "dashed", color = "cornflowerblue") +
  geom_vline(xintercept = -0.1, linetype = "dashed", color = "cornflowerblue") +
  ylim(0, .22) +
  theme_bw()

# More strictly, let's just look at the hypotheses within the 95% CI

eval.rope.mc.cu5 %>%
  filter(grid > -1.226) %>% # This number is the lower bound of the 95% CI: ci(mc.cu4.cu5.post$sample, method = "HDI")$CI_low
  filter(grid < 0.368) %>%  # This number is the upper bound of the 95% CI: ci(mc.cu4.cu5.post$sample, method = "HDI")$CI_high
  pull(rope.prob) %>%
  quantile(probs = seq(.48, .5, length.out = 5))

## Analysis for Study 3 (online data) ----

s3.3s.df <- read.csv("cu7.3s.csv") %>% # Load in Scarlett's online 3 year old data
  as_tibble() %>%
  filter(exclude != "yes") %>% 
  mutate(id = as.numeric(id) + 2000)

# Calculate mean and SD for age
s3.3s.df %>%
  select(birthdate, testdate) %>%
  separate(birthdate, into = c("birth.month", "birth.day", "birth.year")) %>%
  separate(testdate, into = c("test.month", "test.day", "test.year")) %>%
  mutate(birth.year = paste0(20, birth.year),
         test.year = paste0(20, test.year),
         dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.duration(approx.age),
         age.years = age.seconds / (60 * 60 * 24 * 365.25)) %>%
  summarize(mean = mean(age.years),
            sd = sd(age.years))

# How many f / m
table(s3.3s.df$gender)

# Put data in long format
s3.3s.long <- s3.3s.df %>%
  select(id, pretest1:pretest6) %>%
  pivot_longer(-id, names_to = "trial", values_to = "answer") %>%
  filter(!is.na(answer)) %>%
  mutate(trial.type = "3 options",
         cups.or.channels = "slides",
         location = "online",
         response = ifelse(answer == "certain", 1, 0)) %>%
  select(id, cups.or.channels, response)

# Loda Mody & Carey data 
mc.data <- read.csv("modycareydata.csv")

# Put in long format, select the desired trials by age and trial type
mc.3s.choices <- mc.data %>%
  filter(age.group == 3) %>%
  select(id = ID, training1:training3) %>%
  pivot_longer(-id, names_to = "trial", values_to = "response") %>%
  mutate(id = id + 1000,
         cups.or.channels = "cups") %>%
  select(id, cups.or.channels, response)

# Combine the two data frames into one
s3.long <- rbind(mc.3s.choices, s3.3s.long) %>%
  mutate(cups.or.channels = factor(cups.or.channels, levels = c("cups", "slides")))

# Fit a model
# set.seed(3489)
# s3.mod <- stan_glmer(response ~ cups.or.channels + (1|id), family = "binomial", data = s3.long)
# saveRDS(s3.mod, "s3.mod.RDS")
s3.mod <- readRDS("s3.mod.RDS")

s3.means <- emmeans(s3.mod, specs = "cups.or.channels", type = "response")
pairs(emmeans(s3.mod, specs = "cups.or.channels"))

s3.post <- insight::get_parameters(s3.mod) %>%
  as_tibble()

s3.median <- s3.post %>% summarize(median = median(cups.or.channelsslides))
s3.post %>%
  summarize(ci = ci(cups.or.channelsslides, method = "HDI"),
            rope = rope(cups.or.channelsslides, range = c(-0.1, 0.1)))

ggplot(s3.post, aes(x = cups.or.channelsslides)) +
  geom_density()

# Now, let's have a look at all hypotheses in the range of the hypothesis space, and see how 0 stacks up against them.

# This function takes the posterior and calculates the proportion that is between lb and ub (i.e., lower bound and upper bound)
rope.calculator.s3 <- function(lb, ub){
  length(s3.post$cups.or.channelsslides[s3.post$cups.or.channelsslides > lb & s3.post$cups.or.channelsslides < ub]) / nrow(s3.post)
}

# Set up a grid of hypotheses, so that we can evaluate the probability of each hypothesis
eval.rope.s3 <- tibble(grid = seq(min(s3.post$cups.or.channelsslides) %>% round(3), max(s3.post$cups.or.channelsslides) %>% round(3), by = .001), 
                           lower.rope.bound = grid - .1,
                           upper.rope.bound = grid + .1,
                           rope.prob = NA) 

for(i in 1:nrow(eval.rope.s3)){
  eval.rope.s3$rope.prob[i] <- rope.calculator.s3(eval.rope.s3$lower.rope.bound[i], eval.rope.s3$upper.rope.bound[i])
}

eval.rope.s3[eval.rope.s3$grid == 0,] # The probability of the interval around 0 is .236. 
eval.rope.s3[eval.rope.s3$grid == 0.003,] # Here's the probability of the ROPE around the median
quantile(eval.rope.s3$rope.prob, probs = seq(.96, .97, length.out = 5)) # Of all hypotheses, it's better than 96% of hypotheses

ggplot(eval.rope.s3, aes(x = grid, y = rope.prob)) +
  geom_bar(stat = "identity", color = "lightgrey", fill = "lightgrey") +
  geom_vline(xintercept = 0.1, linetype = "dashed", color = "cornflowerblue") +
  geom_vline(xintercept = -0.1, linetype = "dashed", color = "cornflowerblue") +
  ylim(0, .25) +
  theme_bw()

# Find the rank out of all hypotheses in the 95% CI
eval.rope.s3 %>%
  filter(grid > -0.6404303) %>% # This number is the lower bound of the 95% CI: ci(s3.post$cups.or.channelsslides, method = "HDI")$CI_low
  filter(grid < 0.7103767) %>%  # This number is the upper bound of the 95% CI: ci(s3.post$cups.or.channelsslides, method = "HDI")$CI_high
  pull(rope.prob) %>%
  quantile(probs = seq(.92, .93, length.out = 5)) # Top 8% of hypotheses---I'm surprised it doesn't rank better actually

eval.rope.s3 %>%
  filter(rope.prob > .236) %>%
  print(n = 100) # OK, seems legit

# 3-year-old data, first 6 trials only, mean: 60.7%, all trials, 61%

# I don't want to include online 4-year-old data
# Why? They're complicated by experimenter effects; 
# the argument is clear but I don't have space to defend it properly
# 4-year-old data:
# Scarlett's data 50 kids with the first 6 trials only, mean: 65.7%; all trials: 60.5%
# Eimantas' data 26 kids with the first 6 trials only, 50%, all trials, 57% 
# It would be weird to have 74 4-year-olds and 24 3-year-olds

# let's model the data


## Analysis of existing data ----

mc.g.df <- mc.cu.df %>%
  tibble() %>%
  filter(cups.or.channels == "cups" & (age.group == 2 | age.group == 3)) %>%
  mutate(experiment = "mc") %>%
  select(id, age.group, experiment, trial.type, response) %>%
  rbind(g.df %>%
          mutate(experiment = "gr.et.al",
                 age.months = as.numeric(age.months),
                 age.group = factor(ifelse(age.months < 36, 2, 3))) %>%
          select(id, age.group, experiment, trial.type, response)) %>%
  mutate(response = as.numeric(response)) 

mc.g.df$age.group <- droplevels(mc.g.df$age.group)


# Now I'm going to combine all of the data into one dataframe and fit a proper model over this so that I can get pairwise comparisons.
# This is somewhat redundant with the above, but I don't want to delete the code above because I don't know if I would be deleting 
# objects that I would need later on.

full.rep.data <- mc.cu.df %>%
  as_tibble() %>%
  mutate(experiment = ifelse(cups.or.channels == "cups", "Mody & Carey", "Study 1")) %>%
  select(id, age.group, experiment, trial.type, response) %>%
  rbind(cu5.catching.long %>%
          mutate(experiment = "Study 2",
                 age.group = 4) %>%
          select(id, age.group, experiment, trial.type, response)) %>%
  rbind(s3.long %>%
          filter(cups.or.channels == "slides") %>%
          mutate(age.group = 3, 
                 experiment = "Study 3",
                 trial.type = "3options") %>%
          select(id, age.group, experiment, trial.type, response)) %>%
  rbind(mc.g.df %>%
          filter(experiment != "mc"))

# set.seed(48)
# full.mod <- stan_glmer(response ~ age.group * trial.type * experiment + (1|id), family = "binomial", data = full.rep.data)
# saveRDS(full.mod, "full.mod.RDS")
full.mod <- readRDS("full.mod.RDS")

full.mod.means.response <- emmeans(full.mod, specs = "experiment", by = c("age.group", "trial.type"), type = "response") %>%
  as_tibble() %>%
  filter(!is.na(prob)) %>%
  mutate(age.group = factor(age.group, levels = c("2", "3", "4"), labels = c("2.5 yrs", "3 yrs", "4 yrs")),
         trial.type = factor(trial.type, labels = c("3 options", "4 options")),
         experiment = factor(experiment, levels = c("Mody & Carey", "gr.et.al", "Study 1", "Study 2", "Study 3")))

full.mod.pairs <- emmeans(full.mod, specs = "experiment", by = c("age.group", "trial.type")) %>%
  pairs(reverse = TRUE) %>%
  as_tibble() %>%
  filter(!is.na(estimate))

emmeans(full.mod, specs = "trial.type", by = c("age.group", "experiment")) %>%
  pairs(reverse = TRUE) %>%
  as_tibble() %>%
  filter(!is.na(estimate))

full_mod_estimate_means <- full.mod.means.response %>%
  group_by(age.group, trial.type) %>%
  summarize(group_means = mean(prob))

ggplot(full.mod.means.response, aes(x = 1, y = prob, group = experiment)) +
  facet_wrap(~age.group + trial.type, nrow = 1) +
  geom_hline(data = full_mod_estimate_means, aes(yintercept = group_means), linetype = "dashed") +
  geom_errorbar(aes(ymin = lower.HPD, ymax = upper.HPD, color = experiment), position = position_dodge(.35), width = .01) +
  geom_point(aes(color = experiment, shape = experiment), size = 3, fill = "white", position = position_dodge(.35)) +
  scale_shape_manual(values = c(4, 3, 21:23)) +
  ylim(0,1) +  
  theme_bw() +
  theme(text = element_text(size=18), axis.ticks.y=element_blank(), 
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), plot.background = element_blank(), 
        panel.background = element_blank(), legend.title = element_blank(), strip.background = element_blank(),
        axis.ticks.x = element_blank(), axis.text.x = element_blank()) +
  labs(y = "Prob. of choosing target, 95% CI", x = "") 

## Secondary analysis of Mody & Carey (2016) ----

# I need to reload the data to get age as a continuous variable
mc.df.2 <- (read.csv("modycareydata.csv")) # Load mody and carey data

mc.df.long.2 <- mc.df.2 %>%
  select(ID, age..days., age.group, training1:test4) %>%
  pivot_longer(cols = training1:test4, names_to = "trial", values_to = "response") %>%
  mutate(id = ID, age.days = age..days.) %>%
  select(-ID, -age..days.)

mc.df.long.2 <- na.omit(mc.df.long.2)

mc.df.long.2$trial.type <- ifelse(grepl("test", mc.df.long.2$trial), "test", "training")
mc.props.by.kid <- mc.df.long.2 %>%
  group_by(id, age.group, age.days, trial.type) %>%
  summarise(prop.correct = mean(response)) %>%
  mutate(age.group = factor(age.group))
  
# t-tests by age group
# training trials (3 cups)
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "training" & mc.props.by.kid$age.group == 2], mu = 1/3)
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "training" & mc.props.by.kid$age.group == 3], mu = 1/3)
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "training" & mc.props.by.kid$age.group == 4], mu = 1/3)
# test trials
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "test" & mc.props.by.kid$age.group == 2], mu = 1/3)
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "test" & mc.props.by.kid$age.group == 3], mu = 1/3)
t.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "test" & mc.props.by.kid$age.group == 4], mu = 1/3)
# Overall analysis
# 2-way anova
mc.aov <- aov(prop.correct ~ trial.type * age.group, data = mc.props.by.kid)
summary(mc.aov)

# Check assumptions for 2-way ANOVA: Variance Homogenaeity
boxplot(mc.props.by.kid$prop.correct ~ mc.props.by.kid$trial.type * mc.props.by.kid$age.group, main = "prop.correct by trial type × age")
# Boxplots sure look awfully different from one another, looks like a violation?
leveneTest(mc.props.by.kid$prop.correct, interaction(mc.props.by.kid$trial.type, mc.props.by.kid$age.group))
# Levene test comes out nonsignificant, so no violation? But the underlying distribution is binomial, not normal

# Prepare plot
mc.aov.plot <- as.data.frame(allEffects(mc.aov))[[1]]
mc.aov.plot$trial.type <- factor(mc.aov.plot$trial.type, levels = c("training", "test"), labels = c("3opts", "4opts"))
mc.aov.plot$age.group <- factor(mc.aov.plot$age.group, labels = c("2.5yrs", "3yrs", "4yrs", "5yrs"))

ggplot(data = mc.aov.plot) +
  ggtitle("Mody and Carey 2016") +
  facet_wrap(~age.group, nrow = 4) +
  geom_errorbar(mapping = aes(x = trial.type, ymin = lower, ymax = upper), width = .05) +
  geom_point(mapping = aes(x = trial.type, y = fit),
             stat = "identity", alpha = 1, size = 2, shape = 22, fill = "white") +
  ylim(0,1) +
  geom_hline(yintercept = 1/3, linetype = "dashed") +
  theme_bw() +
  theme(text = element_text(size=16, family = "serif"), axis.ticks.y=element_blank(),
        panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), plot.background = element_blank(), 
        panel.background = element_blank(), legend.title = element_blank(), strip.background = element_blank()) +
  ylab("Proportion correct, 95% CIs") +
  xlab("") 

# Age as a continuous variable
# Did they do correlations first?

mc.props.by.kid <- mc.props.by.kid %>%
  mutate(age.days = substr(age.days, 1, nchar(age.days)-1) %>% as.numeric())

cor.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "training"], mc.props.by.kid$age.days[mc.props.by.kid$trial.type == "training"])
cor.test(mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "test"], mc.props.by.kid$age.days[mc.props.by.kid$trial.type == "test"])
# Then for the partial correlation
pcor.mat <- cbind(mc.props.by.kid$age.days[mc.props.by.kid$trial.type == "training"],
                  mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "training"],
                  mc.props.by.kid$prop.correct[mc.props.by.kid$trial.type == "test"])
colnames(pcor.mat) <- c("age.days", "training", "test")

pcor(pcor.mat) # Test the partial correlations. Note that age predicts training trial performance, controlling for test trial performance
# "There was age-related improvement on test trials, over and above that seen on training trials." But there's also
# there was age-related improvement on training trials, over and above that seen on test trials?

# Followup model treating age as a continuous variable
mc.df.long.3 <- mc.df.long.2 %>%
  select(-trial.type) %>%
  pivot_wider(names_from = "trial", values_from = "response") %>%
  mutate(age.standardized = substr(age.days, 1, nchar(age.days)-1) %>% as.numeric %>% scale()) %>%
  rename(age.std = 11) %>%
  pivot_longer(training1:test4, names_to = "trial", values_to = "response") %>%
  mutate(trial.type = ifelse(grepl("training", trial), "training", "test"))
# 
# set.seed(24)
# mc.age.b.mem <- stan_glmer(response ~ trial.type * age.std + (1|id), data = mc.df.long.3, family = binomial, iter = 3500, warmup = 1000)
# saveRDS(mc.age.b.mem, "mc.age.b.mem.RDS")
mc.age.b.mem <- readRDS("mc.age.b.mem.RDS")
summary(mc.age.b.mem)

emtrends(mc.age.b.mem, specs = "trial.type", var = "age.std") %>% pairs()

mc.age.draws <- mc.age.b.mem %>% 
  insight::get_parameters() %>%
  as_tibble() %>%
  rename(interaction = 4)

median(mc.age.draws$interaction)
hdi(mc.age.draws$interaction)

ggplot(mc.age.draws, aes(x = interaction)) +
  geom_density()


## Analysis of age effects ----

# Current 2.5-year-old data
age.data.2.5s <- read.csv("cu3.csv") %>%
  as_tibble() %>%
  filter(exclude != "y") %>%
  separate(birthdate, into = c("birth.month", "birth.day", "birth.year"), sep = "\\/") %>%
  separate(testdate, into = c("test.month", "test.day", "test.year")) %>%
  mutate(birth.year = paste0(20, birth.year),
         birth.month = paste0(0, birth.month),
         birth.day = paste0(0, birth.day),
         birth.day = substr(birth.day, nchar(birth.day) - 1, nchar(birth.day)),
         test.year = paste0(20, test.year),
         test.month = paste0(0, test.month) %>% substr(nchar(test.month), nchar(test.month) + 1),
         test.day = paste0(0, test.day) %>% substr(nchar(test.day), nchar(test.day) + 1),
         dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.duration(approx.age))

age.data.2.5s.long <- age.data.2.5s %>%
  select(id, age.seconds, t3.1:test6) %>%
  mutate(age.seconds = as.numeric(age.seconds)) %>%
  pivot_longer(t3.1:test6, names_to = "trial", values_to = "answer") %>%
  mutate(response = ifelse(answer == "certain", 1, 0),
         trial.type = ifelse(nchar(trial) < 5, "3 options", "4 options"),
         cups.or.channels = "slides",
         blocker = "with.blocker",
         study = "cu.2.5s") %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)

# Study 1 3- and 4-year-old data
age.data.3s.4s <- read.csv("cu4.csv") %>%
  as_tibble() %>%
  filter(exclude == "n") %>%
  separate(birthdate, into = c("birth.month", "birth.day", "birth.year"), sep = "\\/", remove = FALSE) %>%
  separate(testdate, into = c("test.month", "test.day", "test.year"), remove = FALSE) %>%
  mutate(birth.year = paste0(20, birth.year),
         birth.month = paste0(0, birth.month) %>% substr(nchar(birth.month), nchar(birth.month) + 1),
         birth.day = paste0(0, birth.day),
         birth.day = substr(birth.day, nchar(birth.day) - 1, nchar(birth.day)),
         test.year = paste0(20, test.year),
         test.month = paste0(0, test.month) %>% substr(nchar(test.month), nchar(test.month) + 1),
         test.day = paste0(0, test.day) %>% substr(nchar(test.day), nchar(test.day) + 1),
         dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.duration(approx.age))

age.data.3s.4s.long <- age.data.3s.4s %>%
  select(id, age.seconds, t3.1:test6) %>%
  mutate(age.seconds = as.numeric(age.seconds)) %>%
  pivot_longer(t3.1:test6, names_to = "trial", values_to = "answer") %>%
  mutate(response = ifelse(answer == "certain", 1, 0),
         trial.type = ifelse(nchar(trial) < 5, "3 options", "4 options"),
         cups.or.channels = "slides",
         blocker = "with.blocker",
         study = "cu.3s4s") %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)

# Mody & Carey data
mody.carey.age.data <- read.csv("modycareydata.csv") %>%
  as_tibble() %>%
  separate(DOB, into = c("birth.year", "birth.month", "birth.day"), remove = FALSE) %>%
  separate(test.date, into = c("test.year", "test.month", "test.day"), remove = FALSE) %>%
  mutate(dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.duration(approx.age),
         age.years = age.seconds / 3.154e+7,
         id = ID + 1000)

mody.carey.age.data.long <- mody.carey.age.data %>%
  select(id, age.seconds, training1:test4) %>%
  pivot_longer(training1:test4, names_to = "trial", values_to = "response") %>%
  filter(!is.na(response)) %>%
  mutate(trial.type = ifelse(nchar(trial) > 5, "3 options", "4 options"),
         age.seconds = as.numeric(age.seconds),
         cups.or.channels = "cups",
         blocker = "no.blocker",
         study = "mc") %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)

# Study 2 data
cu5.df <- read.csv("cu5.csv") %>% # Load in Eimantas' in-lab 4 year old data
  filter(exclude != "y") %>%
  mutate(age.group = 4) 
cu5.df[cu5.df == "certain "] <- "certain" # Fix a typo

cu5.age.df <- cu5.df %>%   
  separate(birthdate, into = c("birth.month", "birth.day", "birth.year")) %>%
  separate(testdate, into = c("test.month", "test.day", "test.year")) %>%
  mutate(birth.year = paste0(20, birth.year),
         birth.month = paste0(0, birth.month) %>% substr(nchar(birth.month), nchar(birth.month) + 1),
         birth.day = paste0(0, birth.day),
         birth.day = substr(birth.day, nchar(birth.day) - 1, nchar(birth.day)),
         test.year = paste0(20, test.year),
         test.month = paste0(0, test.month) %>% substr(nchar(test.month), nchar(test.month) + 1),
         test.day = paste0(0, test.day) %>% substr(nchar(test.day), nchar(test.day) + 1),
         dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.duration(approx.age))

cu5.age.long <- cu5.age.df %>%
  select(id, age.seconds, pretest1:pretest6) %>%
  pivot_longer(pretest1:pretest6, names_to = "trial", values_to = "answer") %>%
  filter(!is.na(answer)) %>%
  mutate(trial.type = "3 options",
         age.seconds = as.numeric(age.seconds),
         cups.or.channels = "slides",
         blocker = "no.blocker",
         response = ifelse(answer == "certain", 1, 0),
         study = "cu.study2.4s") %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)

# Study 3 data

s3.age.df <- s3.3s.df %>%
  select(id, birthdate, testdate, pretest1:pretest6) %>%
  separate(birthdate, into = c("birth.month", "birth.day", "birth.year")) %>%
  separate(testdate, into = c("test.month", "test.day", "test.year")) %>%
  mutate(birth.year = paste0(20, birth.year),
         test.year = paste0(20, test.year),
         dob = as.Date(paste(birth.year, birth.month, birth.day, sep = "-")),
         dot = as.Date(paste(test.year, test.month, test.day, sep = "-")),
         approx.age = interval(dob, dot),
         age.seconds = as.numeric(as.duration(approx.age))) %>%
  pivot_longer(pretest1:pretest6, names_to = "trial", values_to = "answer") %>%
  mutate(response = ifelse(answer == "certain", 1, 0),
         study = "Study 3",
         trial.type = "3 options",
         cups.or.channels = "slides",
         blocker = "no.blocker") %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)


# Grigoroglou et al. data

g.age.df <- g.df %>%
  mutate(study = "Grigoroglou et al",
         trial.type = factor(trial.type, levels = c("3options", "4options"), labels = c("3 options", "4 options")),
         age.seconds = as.numeric(age.months) * 2628288,
         cups.or.channels = "cups",
         blocker = "no.blocker") %>%
  filter(age.months >= 30) %>%
  select(id, study, age.seconds, trial.type, cups.or.channels, blocker, response)

age.df <- rbind(age.data.2.5s.long, age.data.3s.4s.long, mody.carey.age.data.long, cu5.age.long, s3.age.df, g.age.df) %>%
  mutate(age.years = age.seconds / 31557600,
         response = as.numeric(response)) %>%
  filter(age.years < 5, !is.na(response))

age.df %>%
  group_by(trial.type, cups.or.channels) %>%
  summarize(mean(response))
# Model performance by age, format (cups/slides), and trial type

# This model has all kids younger than 5
# set.seed(804)
# age.mod <- stan_glmer(response ~ age.years * trial.type * cups.or.channels + (1|id), family = "binomial", data = age.df, iter = 3500, warmup = 1000)
# saveRDS(age.mod, "age.mod.RDS")
age.mod <- readRDS("age.mod.RDS")

# This model has all kids except the problematic 4s on 3-slides
# set.seed(804)
# age.mod.2 <- stan_glmer(response ~ age.years * trial.type * cups.or.channels + (1|id), family = "binomial", 
#   data = age.df %>% filter(!(cups.or.channels == "slides" & age.years > 4 & blocker == "with.blocker" & trial.type == "3 options")), 
#   iter = 3500, warmup = 1000)
# saveRDS(age.mod.2, "age.mod.2.RDS")
age.mod.2 <- readRDS("age.mod.2.RDS")

# This model has no slides/cups as a predictor
# set.seed(804)
# age.mod.3 <- stan_glmer(response ~ age.years * trial.type + (1|id), family = "binomial", 
#                         data = age.df, 
#                         iter = 3500, warmup = 1000)
# saveRDS(age.mod.3, "age.mod.3.RDS")
age.mod.3 <- readRDS("age.mod.3.RDS")

# This model has no 4-year-olds at all
# set.seed(804)
# age.mod.4 <- stan_glmer(response ~ age.years * trial.type * cups.or.channels + (1|id), family = "binomial", 
#                         data = age.df %>% filter(age.years < 4), 
#                         iter = 3500, warmup = 1000)
# saveRDS(age.mod.4, "age.mod.4.RDS")
age.mod.4 <- readRDS("age.mod.4.RDS")

age.trends <- emtrends(age.mod, specs = ~ trial.type + cups.or.channels, var = "age.years", type = "response")
age.trends.plot <- as.data.frame(age.trends)

# age.plot <- emmip(age.mod, trial.type + cups.or.channels ~ age.years, at = list(age.years = seq(2.5, 5, by = .05)), type = "response", plotit = FALSE)
# age.plot.2 <- emmip(age.mod.2, trial.type + cups.or.channels ~ age.years, at = list(age.years = seq(2.5, 5, by = .05)), type = "response", plotit = FALSE)
# age.plot.3 <- emmip(age.mod.3, trial.type ~ age.years, at = list(age.years = seq(2.5, 5, by = .05)), type = "response", plotit = FALSE)
# age.plot.4 <- emmip(age.mod.4, trial.type + cups.or.channels ~ age.years, at = list(age.years = seq(2.5, 4, by = .05)), type = "response", plotit = FALSE)

# Participant.props
props.by.age <- age.df %>%
  group_by(id, age.years, trial.type, cups.or.channels) %>%
  summarize(prop.correct = mean(as.numeric(response))) %>%
  filter(age.years >= 2.5)
  

# Let's just quickly have a look at the descriptive data

ggplot(props.by.age, aes(x = age.years, y = prop.correct)) +
  facet_wrap(~ cups.or.channels + trial.type) +
  geom_point(shape = 21) +
  geom_smooth() +
  ylim(0,1) +
  theme_bw() +
  labs(x = "Age", y = "Proportion target choice", title = "All data")

# Descriptives without problematic 4-year-olds
ggplot(age.df %>% 
         filter(!(cups.or.channels == "slides" & age.years > 4 & blocker == "with.blocker" & trial.type == "3 options")) %>%
         group_by(id, trial.type, cups.or.channels, age.years) %>% 
         summarize(prop.correct = mean(response)), 
       aes(x = age.years, y = prop.correct)) +
  facet_wrap(~ cups.or.channels + trial.type) +
  geom_point(shape = 21) +
  geom_smooth() +
  ylim(0,1) +
  theme_bw() +
  labs(x = "Age", y = "Proportion target choice", title = "Excluding problematic 4-year-olds")


# quad.age.mod <- stan_glmer(response ~ (age.years + I(age.years^2)) * trial.type * cups.or.channels + (1|id), family = "binomial", data = age.df)
# saveRDS(quad.age.mod, "quad.age.mod.RDS")
quad.age.mod <- readRDS("quad.age.mod.RDS")

summary(quad.age.mod)

## Guess we're going to have to treat age as categorical

# All data
# set.seed(804)
# age.cat.mod <- stan_glmer(response ~ age.group * trial.type * cups.or.channels + (1|id), family = "binomial",
#                           data = age.df %>%
#                             mutate(age.group = factor(substr(age.years, 1, 1))))
# saveRDS(age.cat.mod, "age.cat.mod.RDS")
age.cat.mod <- readRDS("age.cat.mod.RDS")
# Excluding problematic 4-year-olds
# set.seed(804)
# age.cat.mod.2 <- stan_glmer(response ~ age.group * trial.type * cups.or.channels + (1|id), family = "binomial",
#                           data = age.df %>%
#                             mutate(age.group = factor(substr(age.years, 1, 1))) %>%
#                             filter(!(age.group == 4 & blocker == "with.blocker" & trial.type == "3 options")))
# saveRDS(age.cat.mod.2, "age.cat.mod.2.RDS")
age.cat.mod.2 <- readRDS("age.cat.mod.2.RDS")

# A plot to look at the data: First, all the data
age.cat.points <- pairs(emmeans(age.cat.mod, specs = "age.group", by = c("trial.type", "cups.or.channels")), reverse = TRUE)
age.cat.post <- gather_emmeans_draws(age.cat.points)
# Then the data without the problematic 4s
age.cat.points.2 <- pairs(emmeans(age.cat.mod.2, specs = "age.group", by = c("trial.type", "cups.or.channels")), reverse = TRUE)
age.cat.post.2 <- gather_emmeans_draws(age.cat.points.2)

ggplot(age.cat.post %>% filter(contrast != "4 - 2"), aes(x = .value)) +
  facet_wrap(~cups.or.channels + trial.type + contrast, ncol = 2) +
  geom_density(fill = "lightgrey") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) +
  labs(x = "Effect of age group", y = "Density")

ggplot(age.cat.post.2 %>% filter(contrast != "4 - 2"), aes(x = .value)) +
  facet_wrap(~cups.or.channels + trial.type + contrast, ncol = 2) +
  geom_density(fill = "lightgrey") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) +
  labs(x = "Effect of age group", y = "Density")



## Further analysis: effect of trial number, using only non-artifacted, in-lab data ----

no.artifact.inlab.df <- cu5.catching.long %>%
  mutate(age.group = 4,
         cups.or.channels = "slides",
         trial.number = as.numeric(substr(trial, nchar(trial), nchar(trial)))) %>%
  select(id, age.group, cups.or.channels, trial.type, trial.number, response) %>%
  rbind(mc.cu.df %>%
          select(id, age.group, cups.or.channels, trial.type, trial.number, response) %>%
          filter(!(age.group == 4 & cups.or.channels == "slides" & trial.type == "3options")))

# set.seed(7492)
# no.artifact.inlab.trial.num.mem <- stan_glmer(response ~ cups.or.channels * trial.type * age.group * trial.number + (1|id),
#                                   family = "binomial", data = no.artifact.inlab.df, iter = 3500, warmup = 1000)
# saveRDS(no.artifact.inlab.trial.num.mem, "no.artifact.inlab.trial.num.mem.RDS")
no.artifact.inlab.trial.num.mem <- readRDS("no.artifact.inlab.trial.num.mem.RDS")

no.artifact.inlab.trial.num.means <- emtrends(no.artifact.inlab.trial.num.mem, specs = ~ cups.or.channels + trial.type + age.group, var = "trial.number", type = "response")
no.artifact.inlab.trial.num.post <- gather_emmeans_draws(no.artifact.inlab.trial.num.means)

no.artifact.inlab.trial.num.plot <- emmip(no.artifact.inlab.trial.num.mem, cups.or.channels + trial.type + age.group ~ trial.number, at = list(trial.number = seq(1, 6, by = .5)), CIs = TRUE, plotit= FALSE, type = "response")

no.artifact.trial.number.props <- no.artifact.inlab.df %>%
  filter(cups.or.channels == "slides") %>%
  group_by(age.group, trial.type, trial.number) %>%
  summarize(prop.correct = mean(response))

no.artifact.cis <- no.artifact.inlab.trial.num.post %>%
  filter(cups.or.channels == "slides") %>%
  group_by(age.group, trial.type) %>%
  summarize(ci = ci(.value, method = "HDI")) %>%
  unnest(ci)

ggplot(no.artifact.inlab.trial.num.plot %>% filter(cups.or.channels != "cups"), aes(x = trial.number, y = yvar)) +
  facet_wrap(~age.group + trial.type, nrow = 3) +
  geom_smooth() +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = .25, color = "white") +
  # geom_point(data = no.artifact.trial.number.props, aes(y = prop.correct)) +
  ylim(0,1) +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

ggplot(no.artifact.inlab.trial.num.post %>% filter(cups.or.channels == "slides"), aes(x = .value)) +
  facet_wrap(~age.group + trial.type, nrow = 3) +
  geom_density(fill = "lightgrey") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_segment(data = no.artifact.cis, aes(x = CI_low, y = 0, xend = CI_low, yend = 1.5), color = "darkred") +
  geom_segment(data = no.artifact.cis, aes(x = CI_high, y = 0, xend = CI_high, yend = 1.5), color = "darkred") +
  theme_bw()

# What if I include all of the unbiased slides data?

no.artifact.slides.df <- cu5.catching.long %>%
  mutate(age.group = 4,
         cups.or.channels = "slides",
         trial.number = as.numeric(substr(trial, nchar(trial), nchar(trial)))) %>%
  select(id, age.group, cups.or.channels, trial.type, trial.number, response) %>%
  rbind(mc.cu.df %>%
          select(id, age.group, cups.or.channels, trial.type, trial.number, response) %>%
          filter(!(age.group == 4 & cups.or.channels == "slides" & trial.type == "3options"), cups.or.channels != "cups")) %>%
  rbind(s3.3s.df %>%
          select(id, pretest1:pretest6) %>%
          pivot_longer(-id, names_to = "trial", values_to = "answer") %>%
          filter(!is.na(answer)) %>%
          mutate(age.group = "3",
                 trial.type = "3options",
                 cups.or.channels = "slides",
                 location = "online",
                 response = ifelse(answer == "certain", 1, 0),
                 trial.number = substr(trial, nchar(trial), nchar(trial)) %>% as.numeric()) %>%
          select(id, age.group, cups.or.channels, trial.type, trial.number, response))

no.artifact.slides.props <- no.artifact.slides.df %>%
  group_by(age.group, trial.type, trial.number) %>%
  summarize(prop.correct = mean(response)) %>%
  as_tibble()

# set.seed(7492)
# no.artifact.slides.trial.num.mem <- stan_glmer(response ~  trial.type * age.group * trial.number + (1|id),
#                                   family = "binomial", data = no.artifact.slides.df, iter = 3500, warmup = 1000)
# saveRDS(no.artifact.slides.trial.num.mem, "no.artifact.slides.trial.num.mem.RDS")
no.artifact.slides.trial.num.mem <- readRDS("no.artifact.slides.trial.num.mem.RDS")

no.artifact.slides.trial.num.trends <- emtrends(no.artifact.slides.trial.num.mem, specs = ~ trial.type + age.group, var = "trial.number", type = "response")
no.artifact.slides.trial.num.post <- gather_emmeans_draws(no.artifact.slides.trial.num.trends)
no.artifact.slides.trial.num.plot <- emmip(no.artifact.slides.trial.num.mem, trial.type + age.group ~ trial.number, at = list(trial.number = seq(1, 6, by = .5)), CIs = TRUE, plotit= FALSE, type = "response") %>%
  as_tibble()


ggplot(no.artifact.slides.trial.num.plot %>%
         mutate(age.group = paste(age.group, "yrs"),
                trial.type = ifelse(trial.type == "3options", "3 options", "4 options")), aes(x = trial.number, y = yvar)) +
  facet_wrap(~age.group + trial.type, nrow = 3) +
  geom_smooth(color = "darkgrey", width = .75) +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = .5) +
  geom_point(data = no.artifact.slides.props %>%
               mutate(age.group = paste(age.group, "yrs"),
                      trial.type = ifelse(trial.type == "3options", "3 options", "4 options")), aes(x = trial.number, y = prop.correct)) +
  ylim(0,1) +
  theme_bw() +
  labs(x = "Trial number", y = "Probability of choosing target")
  
no.artifact.slides.trial.num.post %>%
  group_by(trial.type, age.group) %>%
  summarize(median = median(.value))

no.artifact.slides.trial.num.post %>%
  group_by(age.group, trial.type) %>%
  summarize(ci = ci(.value, method = "HDI")) %>%
  unnest(ci)


## Individual-level (mixture) analysis ----
## Setup, define functions

noisy.prop.prob.calculator <- function(estimate, noise.rate){
  prob = seq(0,1, length.out = 1001)
  prop.prob = (estimate - noise.rate) / (prob - noise.rate)
  return(tibble(prob = prob) %>%
           mutate(prop.probabilistic = ifelse(prop.prob <= 0, 0, ifelse(prop.prob >= 1, 1, prop.prob)),
                  prop.noisy.max = 1-prop.probabilistic,
                  noise.rate = noise.rate))
}


noisy.probability.observed.from.mix.calculator <- function(observed.vector, tib, prior){
  n = sum(observed.vector)
  local.prob <- vector()
  global.prob <- vector()
  numerator <- vector()
  for(j in 1:nrow(tib)){
    for(i in 1:length(observed.vector)){ 
      
      local.prob[i] <- dbinom(nth(observed.vector, i), n, dbinom(i-1, length(observed.vector) - 1, tib$prob[j]) * tib$prop.probabilistic[j] + 
                                dbinom(i - 1, length(observed.vector) - 1, tib$noise.rate[j]) * tib$prop.noisy.max[j])
    }
    global.prob[j] <- prod(local.prob)
    numerator[j] <- global.prob[j] * prior$prob[j]
  }
  
  return(tib %>%
           mutate(likelihood = global.prob,
                  probability.of.observed.dist.under.this.mix = numerator,
                  normalizing.constant = sum(numerator),
                  normalized.probability = numerator / normalizing.constant))
}

# This function calculates 95% HDIs
my.95.hdi <- function(data){
  id.vec <- NA
  sum.vec <- NA
  above.95 <- FALSE
  j <- 1
  while(!above.95){
    sum.vec[j] <- max(data$normalized.probability[!data$prob %in% id.vec])
    id.vec[j] <- data %>% filter(!prob %in% id.vec) %>% filter(normalized.probability == max(normalized.probability)) %>% pull(prob)
    j <- j+1
    above.95 <- ifelse(sum(sum.vec) > .95, TRUE, FALSE)
  }
  c(min(id.vec), max(id.vec))
}


## OK Now the actual tests
# Observed distributions:
slides.dist.3options.2s <- readRDS("slides.dist.3options.2s.RDS")
slides.dist.3options.3s <- readRDS("slides.dist.3options.3s.RDS")
slides.dist.3options.4s <- readRDS("slides.dist.3options.4s.RDS")
slides.dist.4options.2s <- readRDS("slides.dist.4options.2s.RDS")
slides.dist.4options.3s <- readRDS("slides.dist.4options.3s.RDS")
slides.dist.4options.4s <- readRDS("slides.dist.4options.4s.RDS")

prior <- tibble(x = seq(0, 1, length.out = 1001),
                prob = (dnorm(seq(0, 1, length.out = 1001), mean = 1/3, sd = 0.02) + dnorm(seq(0, 1, length.out = 1001), mean = 1/2, sd = 0.02)) / 2)

ggplot(prior, aes(x = x, y = prob)) +
  geom_line() +
  geom_vline(xintercept = 1/3, linetype = "dashed") +
  geom_vline(xintercept = 1/2, linetype = "dashed") +
  theme_bw() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), panel.grid = element_blank()) +
  labs(x = "Hypothesized ground-truth probability among nonmaximizers", y = "Prior probabillity")

area_under_curve(x = prior$x, y = prior$prob)
## Model existing data

# set.seed(824)
# age.all.trial.type.mem <- stan_glmer(response ~ age.group * trial.type + (1|id), family = "binomial",
#                           data = no.artifact.inlab.df %>% filter(cups.or.channels == "slides"),
#                           iter = 3500, warmup = 1000)
# saveRDS(age.all.trial.type.mem, "age.all.trial.type.mem.RDS")
age.all.trial.type.mem <- readRDS("age.all.trial.type.mem.RDS")
# set.seed(824)
# age.all.mem <- stan_glmer(response ~ age.group + (1|id), family = "binomial", 
#                           data = no.artifact.inlab.df %>% filter(cups.or.channels == "slides"), 
#                           iter = 3500, warmup = 1000)
# saveRDS(age.all.mem, "age.all.mem.RDS")
age.all.mem <- readRDS("age.all.mem.RDS")

# loo1 <- loo(age.all.trial.type.mem)
# loo2 <- loo(age.all.mem)
# loo_compare(loo1, loo2) # LOO prefers the simpler model

age.all.means <- emmeans(age.all.mem, specs = "age.group", type = "response")

# > age.all.means
# age.group  prob lower.HPD upper.HPD
# 2         0.439     0.340     0.542
# 3         0.604     0.498     0.705
# 4         0.696     0.612     0.776

# Or maybe just use the descriptive means of the unbiased inlab data:
no.artifact.inlab.df %>% 
  filter(cups.or.channels == "slides") %>%
  group_by(age.group) %>%
  summarize(prop.correct = mean(response))

## Noise rate = 1
## First, 2s
age.noise.1.3opts.2s <- noisy.prop.prob.calculator(estimate = .446, noise.rate = 1) %>%
  mutate(age.group = 2, trial.type = "3options")
age.noise.1.4opts.2s <- noisy.prop.prob.calculator(estimate = .446, noise.rate = 1) %>%
  mutate(age.group = 2, trial.type = "4options")

saved.vals.age.noisy.1.3opts.2s <- noisy.probability.observed.from.mix.calculator(slides.dist.3options.2s, age.noise.1.3opts.2s, prior)
saved.vals.age.noisy.1.4opts.2s <- noisy.probability.observed.from.mix.calculator(slides.dist.4options.2s, age.noise.1.4opts.2s, prior)

## Next, 3s
age.noise.1.3opts.3s <- noisy.prop.prob.calculator(estimate = .594, noise.rate = 1) %>%
  mutate(age.group = 3, trial.type = "3options")
age.noise.1.4opts.3s <- noisy.prop.prob.calculator(estimate = .594, noise.rate = 1) %>%
  mutate(age.group = 3, trial.type = "4options")

saved.vals.age.noisy.1.3opts.3s <- noisy.probability.observed.from.mix.calculator(slides.dist.3options.3s, age.noise.1.3opts.3s, prior)
saved.vals.age.noisy.1.4opts.3s <- noisy.probability.observed.from.mix.calculator(slides.dist.4options.3s, age.noise.1.4opts.3s, prior)

## Last, 4s
age.noise.1.3opts.4s <- noisy.prop.prob.calculator(estimate = .667, noise.rate = 1) %>%
  mutate(age.group = 4, trial.type = "3options")
age.noise.1.4opts.4s <- noisy.prop.prob.calculator(estimate = .667, noise.rate = 1) %>%
  mutate(age.group = 4, trial.type = "4options")

saved.vals.age.noisy.1.3opts.4s <- noisy.probability.observed.from.mix.calculator(slides.dist.3options.4s, age.noise.1.3opts.4s, prior)
saved.vals.age.noisy.1.4opts.4s <- noisy.probability.observed.from.mix.calculator(slides.dist.4options.4s, age.noise.1.4opts.4s, prior)

age.noise.1 <- rbind(saved.vals.age.noisy.1.3opts.2s, saved.vals.age.noisy.1.3opts.3s, saved.vals.age.noisy.1.3opts.4s,
                     saved.vals.age.noisy.1.4opts.2s, saved.vals.age.noisy.1.4opts.3s, saved.vals.age.noisy.1.4opts.4s) 


# Calculate 95% CIs

all.saved.vals.cis <- saved.vals.age.noisy.1.3opts.2s %>%
  mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.2s)[1],
         max.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.2s)[2]) %>%
  rbind(saved.vals.age.noisy.1.3opts.3s %>%
          mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.3s)[1],
                 max.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.3s)[2])) %>%
  rbind(saved.vals.age.noisy.1.3opts.4s %>%
          mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.4s)[1],
                 max.95 = my.95.hdi(saved.vals.age.noisy.1.3opts.4s)[2])) %>%
  rbind(saved.vals.age.noisy.1.4opts.2s %>%
          mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.2s)[1],
                 max.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.2s)[2])) %>%
  rbind(saved.vals.age.noisy.1.4opts.3s %>%
          mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.3s)[1],
                 max.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.3s)[2])) %>%
  rbind(saved.vals.age.noisy.1.4opts.4s %>%
          mutate(min.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.4s)[1],
                 max.95 = my.95.hdi(saved.vals.age.noisy.1.4opts.4s)[2]))

all.saved.vals.summary <- all.saved.vals.cis %>%
  group_by(age.group, trial.type) %>%
  summarize(lcl = mean(min.95),
            ucl = mean(max.95))

# This code calculates the rank of .5 and .33 among the hypotheses in the 95% CI
all.saved.vals.inside.ci <- all.saved.vals.cis %>%
  mutate(too.low = ifelse(min.95 > prob, 1, 0),
         too.high = ifelse(max.95 < prob, 1, 0)) %>%
  filter(too.low == 0 & too.high == 0) %>%
  group_by(age.group, trial.type, noise.rate) %>%
  arrange(desc(normalized.probability), by_group = TRUE) %>%
  mutate(rank = row_number(),
         quantile = rank / max(rank))

all.saved.vals.inside.ci %>% filter(prob == .5)
all.saved.vals.inside.ci %>% filter(prob == .33)

# Visualize the likelihood
ggplot(age.noise.1, aes(x = prob, y = likelihood)) +
  facet_wrap(~age.group + trial.type, ncol = 2, scales = "free_y") +
  geom_line(size = .75) +
  geom_area(fill = "lightgrey") +
  geom_vline(xintercept = .5, linetype = "dashed") +
  geom_vline(xintercept = 1/3, linetype = "dashed") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        strip.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.spacing.y = unit(2, "lines"),
        text = element_text(size = 18)) +
  labs(x = "Likelihood: pr(data | hypothesis)",
       y = "")

# Visualize the posterior
ggplot(age.noise.1, aes(x = prob, y = normalized.probability)) +
  facet_wrap(~age.group + trial.type, ncol = 2) +
  geom_line(size = .75) +
  geom_area(fill = "lightgrey") +
  geom_vline(xintercept = .5, linetype = "dashed") +
  geom_vline(xintercept = 1/3, linetype = "dashed") +
  geom_segment(data = all.saved.vals.summary, aes(x = lcl, xend = lcl, y = 0, yend = .01), size = .5, color = "darkred") +
  geom_segment(data = all.saved.vals.summary, aes(x = ucl, xend = ucl, y = 0, yend = .01), size = .5, color = "darkred") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        strip.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.spacing.y = unit(2, "lines"),
        text = element_text(size = 18)) +
  labs(x = "Ground-truth probability of choosing target in non-maximizers",
       y = "Probability")

## Which 2.5-year-olds performed well on 4-slides? How did they perform on 3-slides and 2-slides?
# Find their ids
perfect.ids.2.5.4slides <- mc.cu.df %>%
  as_tibble() %>%
  filter(age.group == 2 & trial.type == "4options" & cups.or.channels == "slides") %>%
  group_by(id) %>%
  summarize(num.correct = sum(response)) %>%
  filter(num.correct == 6) 

# Find what their prop correct was on 2 slides
cu.2cups.long %>%
  filter(id %in% perfect.ids.2.5.4slides$id) %>%
  group_by(id) %>%
  summarize(prop.correct = mean(response))

# Find what their prop correct was on 3 slides
mc.cu.df %>%
  filter(id %in% perfect.ids.2.5.4slides$id & trial.type == "3options") %>%
  group_by(id) %>%
  summarize(prop.correct = mean(response))



